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Abstract 

The paper addresses the generalization of the half-quadratic minimization method 
for the restoration of images having values in a complete, connected Riemannian 
manifold. We recall the half-quadratic minimization method using the notation of 
the c-transform and adapt the algorithm to our special variational setting. We prove 
the convergence of the method for Hadamard spaces. Extensive numerical examples 
for images with values on spheres, in the rotation group SO(3), and in the manifold 
of positive definite matrices demonstrate the excellent performance of the algorithm. 
In particular, the method with SO(3)-valued data shows promising results for the 
restoration of images obtained from Electron Backscattered Diffraction which are of 
interest in material science. 


1 Introduction 

Many edge-preserving variational methods for the denoising or inpainting of real-valued 
images utilize the following model: let Q := {1,... ,n} x ,m} be the image grid, 

0 ^ V C ^ and M{i)^ {(ii + 1, ^ 2 ), (m, '^2 + 1)} the set of right and upper neighbors of 
pixel i ^ where we suppose mirror boundary conditions. From corrupted image values 
/: V ^ M we want to restore the original image ito: 5 ^ M as a minimizer of one of the 
following energy functionals 


- Uif + XI (1) 

lev ieG 



where A > 0 is a regularization parameter and (f\ M>o ^ M>o. For V = Q this is a typical 
denoising model in the presence of additive Gaussian noise. Otherwise, the model can be 
used for inpainting the missing image values in Q\y. Throughout this paper, we consider 
even functions (^: M ^ M>o. For ip{t) := |t|, the models (1) and (2) are discrete variants 
of the anisotropic and isotropic Rudin-Osher-Fatemi model [42], respectively. Then the 
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regularization term is often referred as anisotropic/isotropic discrete total variation (TV) 
regularization in resemblance to its functional analytic counterpart. 

Remark 1.1. More generally one may eonsider Q as vertiees of a graph with edge set 
£ {{hj) • i ^ ^ •^(^)} some appropriate neighborhoods Af{i). This is for 

example useful in nonloeal means approaehes. Then the regularizing term sums over the 
edge set £. For simplieity we restriet our attention to the speeial neighborhoods J\f{i)^ 
here. 

Starting with the inaugural work [21,22], a huge number of papers have examined the 
so-called half-quadratic minimization methods for solving the above restoration problems 
with various functions (p as well as other optimization problems. We only mention the 
ARTUR algorithm in [16]. Basically, the original problem is reformulated into an aug¬ 
mented one which is quadratic with respect to the image and separable with respect to 
additional auxiliary variables. Then an alternating minimization process is applied whose 
steps allow an efficient computation. Half-quadratic minimization is connected with other 
well-known minimization approaches. We only mention the relation to EM algorithms [14], 
quasi-Newton minimization [4,33] and gradient linearization algorithms [33]. A gradient 
linearization method was used in particular in [53,54] to minimize an approximate total 
variation regularization (2) with (p{t) := ^/W+~^ for 6: <C 1. It was called the “lagged 
diffusivity fixed point iteration” and the authors mention that it amounts to apply the 
multiplicative form of half-quadratic minimization to this (p. Finally, there is a relation to 
iteratively reweighted least squares methods [18,30]. For a convergence analysis of half¬ 
quadratic minimization methods for convex functions p) we refer to [16,34] and for weaker 
convergence results for nonconvex p) to [19]. 

In many applications signals or images having values in a manifold are of interest. 
Circle-valued images appear in interferometric synthetic aperture radar [13,20] and various 
applications involving the phase of Fourier transformed data. Images with values in 
play a role when dealing with 3D directional information [27,29,51] or in the processing 
of color images in the chromaticity-brightness (CB) setting [15]. The motion group and 
the rotation group SO(3) were considered in tracking, (scene) motion analysis [38,40,50] 
and in the analysis of back scatter diffraction data [8]. Finally, images with values in 
the manifold of positive definite matrices appear in DT-MRI [36,45,55,57] and whenever 
covariance matrices are adjusted to image pixels, see, e.g., [50]. 

Recently a TV-like model for circle-valued images was introduced in [46,47]. For 
manifold-valued image restoration such an approach was proposed in [31], where the prob¬ 
lem was reformulated as a multilabel optimization problem which was handled using con¬ 
vex relaxation techniques. Another method suggested in [56] employs cyclic and parallel 
proximal point algorithms and does not require labeling and relaxation techniques. This 
approach was generalized by including second order differences for circle-valued images 
in [9,10], for coupled circle and real-valued images in [11] and for Riemannian manifolds 
in [6]. A restoration method which circumvents the direct work with manifold-valued data 
by embedding the matrix manifold in the appropriate Euclidean space and applying a back 
projection to the manifold was suggested in [41]. Recently, an iteratively reweighted least 
squares method for the restoration of manifold-valued images was suggested in [24]. This 
method can be seen as multiplicative half-quadratic minimization method for the special 
function p){t) := 

In this paper we adopt the idea of half-quadratic minimization for general functions p) 
for the restoration of manifold-valued images. We prefer the notation of the c-transform 
known from optimal transport to recall the basic half-quadratic minimization approach. 
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Then we describe the algorithm for our problems of denoising or inpainting of manifold¬ 
valued images both in the anisotropic and isotropic case. Here we focus on the multiplica¬ 
tive half-quadratic minimization method. Convergence of the algorithm can be shown 
for images having entries in an Hadamard space. The manifold of positive definite ma¬ 
trices is such an Hadamard manifold. We provide several applications of the algorithm 
as the denoising of phase-valued images, the restoration of color images with disturbed 
chromaticity or of 3D directions, and the improvement of images obtained from electron 
backscatter diffraction of a Magnesium sample. 

The outline of the paper is as follows: In Section 2 we propose our variational model and 
show how to handle it by the half-quadratic minimization approach. A convergence proof 
for the algorithm applied on Hadamard manifolds is given in Section 3. Section 4 shows 
various numerical examples. Finally, Appendix A contains the proofs and Appendix B 
lists special quantities necessary for the numerical computations. 

2 Half-Quadratic Minimization 

Let A4 be a complete, connected n-dimensional Riemannian manifold with geodesic dis¬ 
tance d: A4 X A4 ^ M>o. Now we consider manifold-valued images. More precisely, from 
corrupted image values /: V ^ A4 we want to restore the original manifold-valued image 
0 ^ A4 as a minimizer of one of the following energy functionals 

•= + XI (3) 

iev ieg j&j\f(i)+ 

Mu) X d^{ui,Uj)y'] (4) 

iev ieG ^ jGA/'(b+ ^ 

with A > 0 being a regularization parameter and R>o ^ R>o- As before we set 
(p{t) := ^{—t) for t < 0 and consider as a function defined on the whole real axis. For 
ip{t) := |t|, the second sum in the regularizer of Jy is just || {d{ui^ jGA/'(i)+ IL’ ^ ^ 

Then Ji resembles the setting in [56] and is related to the anisotropic ROF functional 
(1) and J 2 gives the approach [31] related to the isotropic case (2). In this paper we 
will consider smooth regularization terms, i.e., even, differentiable functions ip. We will 
compute a minimizer of J^y, v G {1,2}, by half-quadratic minimization methods. 

In the following, we briefly recall the reformulation idea of half-quadratic minimization 
using the concept of the c-transform and apply it to (3) and (4). The c-transform of 
functions defined on metric spaces is used in connection with optimal transport problems 
see, e.g., [52, p. 86f] and seems also to be an appropriate approach here. Given a function 
c : M X M ^ M, the c-transform of a function (^ : M ^ M is defined by 

:= inf{c(t,s) - 

By this definition we see that ip{t)+(p^{s) < s). For c{t^ s) := —st we have (p^ = —(—(^)* 

with the Fenchel transform defined for /i : M ^ M as 

/i*(5) := sup{t5 — h{t)Y 

Recall that h** — h if and only if h is lower semi-continuous (Isc) and convex. In the 


3 



multiplicative and additive half-quadratic methods we use the functions 


c(t, s) := t^s, 


(multiplicative) 

(5) 

c(t,s) := i( 

^ 1 

Va t - ^s] , 

^fa J 

a > 0, (additive) 

(6) 


respectively. The multiplicative setting was introduced in [21] and the additive one in [22]. 
The quadratic cost function c(t, s) in the additive setting is also handled in optimal trans¬ 
port topics. Note that the cost function c(t, s) in the multiplicative case is not bounded 
from below in s. The following proposition is crucial for the half-quadratic reformulation. 

Proposition 2.1. Let (^: M ^ M>o be an even, differentiable function and /et c: M x M ^ 
M be given by (5) and (6)^ respectively. 

i) If the functions 


$(t) := 

:= 


fort>0, 
^+oo for t < 0, 


(multiplicative) 

(additive) 


respectively, are convex, then ip — (p^^, i.e., setting fj{s) := ip^{s) one has 


= in|{c(i,s) - V’(s)}- 


(7) 

( 8 ) 


ii) If in addition to the assumption in i) we have 

lim ^ 0, (multiplicative) (9) 

t^OO t^ 

lim < ^a, (additive) (10) 

t^OO 2 

respectively, and in the multiplicative case also ip'{t) > 0 for t > 0 and (/^"(0-h) := 
limt^o+ exists, then the infimum in (7) and (8) is attained for {t,s) — {t,s{t)) 


s{t) := 


^ for t > 0, 


fort = 0, 

s{t) := at — ip'{t), 


(multiplicative) 


(additive) 


respectively, and for these pairs we have p){t) + f>{s) = c{t, 5). The choice is unique 
except for the multiplicative case and t = 0, where any s larger than ^ is also a 
solution. 


hi) If in the multiplicative case in addition p>\t) > 0 for t > 0 and (^"(0-h) > 0^ then 


Note that in the multiplicative case fj{s) — —00 for 5 < 0, so that we can restrict our 
attention in the infimum in (8) to 5 > 0. Further, the assumption p^\t) > 0, t > 0, is in 
particular fulfilled if p){t) is convex and (^"(0-h) > 0. The proof can be given following for 


4 



example the lines in [16,34]. However, since the assumptions in these papers are slightly 
different and the c-transform notation is not used, we add the proof in the appendix A to 
make the paper self-contained. The functions ip which fulfill the conditions in Proposition 
2.1 and which were used in our numerical test are listed in Table 1. Further examples are 
collected in [16,33,34]. 

In the following we assume that p) fulfills the assumptions of Proposition 2.1. Now the 
idea is to replace p in (3), resp., (4) by the expression in (8) and to consider min,^ Ju{u) = 
rninu min-i; v) with 


Jl{u,v) : 




(12) 


iev 

ieG jGA/'(i)+ 


J2{u,v) : 


^*) + aE( 

^( ( XI d^{ui,Uj)Y ,Vi\ 

, (13) 


iev 

ieg 

V A A J 



where we have used the notation v := in the anisotropic case and v := {vi)ieg 

in the isotropic case. By Proposition 2.1 i), minimizing over u and v gives the same 
solutions for u as just minimizing over u. More precisely we give the following remark. 

Remark 2.2. Let us abbreviate d^j := d{ui^Uj), d-^ := in the anisotropie ease 

1 ’ 

and d^ := ( ^ d‘^{ui^Uj))‘^, d-^ := {di)i^g in the isotropie ease. If u is a minimizer of 

Jjy, z/ = 1,2, then (^u^s{du)) is a minimizer of z/ = 1,2 and eonversely. In partieular, 

if 

u — argmin s{u)) 

u 

holds true, then u is a minimizer of J^. 

Now we can apply an alternating minimization over vi^j G M, resp., G M and u G 


^(^+1) ^ argmin 'z;), (14) 

V 

^(^+1) ^ argmin (15) 

u 

Clearly, under the assumptions of Proposition 2.1 we have 

jpui’‘\vi’^+^^)=jpui’^^). (16) 

We want to work with the differentiable function d^ in the second iteration (15). Since 
the additive reformulation leads to a non-differentiable function d^ + sd, we restrict our 
attention in the following to the multiplicative case. 


Minimization with respect to v. The minimization over v in (14) can be done sepa¬ 
rately for the Vij or Vi. By Proposition 2.1 a minimizer is given by 

:= s(d„(fc)), 

where s is defined as in (11). Note that only for d{pp\uP^ = 0 in the anisotropic case 

1 

and ^ ^ = 0 in the isotropic case a larger value also could be taken 

as a minimizer. 

If p fulfills the assumptions of Proposition 2.1 iii), then € (0,<^"(0+)/2]. 
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Minimization with respect to u. The minimization over u in (15) is equivalent to 
finding the minimize! of 

(«) + A 

iev i£g j&J\f{i)+ 

i^V i^Q jGA/'(i)+ 

respectively. We can apply, e.g., a gradient descent or a Riemann-Newton method, see [1]. 
Both methods are described in the following for our setting and were implemented. 

We need the following notation. Let TxM denote the tangential space of Ad at x G Ad 
and (*, •)a;: TxAi x TxM M the Riemannian metric with induced norm ll-Ha,. Let 
X G Ad, G TxM be the minimal geodesic starting from 7a;,^(0) = x with 7a;,^(0) = C* 
Then the exponential map exp^: TxM. ^ Ad is given by exp^(^ = 7a;,^(l). The inverse 
exponential map denoted by log^, = exp“^: Ad ^ TxM is locally well-defined. For the 
manifolds used in our numerical examples, namely the sphere n G N, the SO(3) and 
the manifold P(r), r G N, of symmetric positive definite r x r matrices, the specific maps 
are given in the Appendix B. Finally, for F: Ad ^ M, let 

gradF(x) G TxM and Hessi;^(x) : TxM TxM 

be the Riemannian gradient and the Hessian ofFiAd^MatxGAd, respectively. For 
(i(*, y) : M ^ M>o one has 

gradcf{x,y) = -21og3.y. 

Considering the image G M, we abbreviate the gradient and the Hessian of a function 
F : M ^ M defined on the product manifold M also by gradF and Hessi^^, resp., since its 
use becomes clear from the context. 

Gradient descent method. The gradient descent method computes, starting with 
ui^) iteratively 

«(’■+!) = exp-M (trgrad v G {1,2}, (17) 

with appropriate step sizes F > 0. The gradient grad j7’^^(/c+i) is given by 

(grad (m)) . = -lv(i) log„. /* - 2A ^ (18) 

(gradJ 2 yfc+i)(’^)). =-lv(i)log„/i-2A(?;f+^^ vf^'^hog^u}j (19) 

jGA/'(i)+ jeU{i)- 

where i G 0 , M{i)~ := {(ii - IA2), (nA2 - 1)} and := M{i)^ 

Riemann-Newton method. Alternatively we can use a Riemann-Newton method to 
compute a minimizer. Finding a descent direction r]r ^ with Newton’s method is 

done for z/ G {1, 2} by solving the system of equations 

(«(’■)) (rjr) = - grad J^^^(k+i) (20) 
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Algorithm 1 Image Restoration by Half-Quadratic Minimization (multiplicative) 

Input: V, corrupted image / GAd^^CM,A,(^ 

Output: Restored image u 

Initialize 

repeat 

k i — k 

vik+l) ^ 

Compute 

uik+i) exp^(R) riR 

by R steps of a gradient descent method (18), resp. (19) 
or by Newton’s approach (20); 
until stopping criterion is reached; 


Then we update, starting with := iteratively 

dr+l) _ 




exp~(.) r]r 
exp~(.) 


if (r?r,grad < 0, 

(— grad J^^^(k+i) ('W*-’^^)) otherwise. 


The whole half-quadratic minimization method for our problem is given in Algorithm 1. 


In [33] it was shown that the multiplicative half-quadratic minimization is equivalent to 
the quasi-Newton descent method, and therefore we expect that its performance is better 
than the simple gradient descent method. Let us comment on this for the manifold-valued 
setting. 


Remark 2.3. (Relation of half-quadratic minimization to gradient descent and (quasi) 
Newton methods) 

We restrict our attention to the case v — 1. Similar considerations can he done for v — 2. 
The gradient of the initial functional Ji in (3) is given for ui ^ Uj by 


(grad Ji(n)). 


-lv(i) log„,/i - A 

jeAf(i) 


-lv(i) /i - 2A 


ip'{d{ui,Uj)) 


2d(ui,Un) 

jeAfii) 


log«, Uj 

log«i Uj 


which by (11) can he rewritten as 

(grad = -lv(i) log„, fi-2\ ^ s{d{ui,Uj)) log„, uj. (21) 

jeN{i) 

Hence a gradient descent algorithm applied to the initial functional Ji coincides with the 
half-quadratic method if we perform only one step in the gradient descent method (17) to 
obtain an update of u. More than one gradient descent step in (17) leads to a linearized 
gradient descent of R. If we perform a Riemann-Newton step (20) to update u we have 
a quasi-Newton method for R. Note that — s{d^(k)) is fixed in the half-quadratic 

update step of u which is not the case in (21). This simplifies the computation of the 
Hessian in the half-quadratic approach. 
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3 Convergence in Hadamard Manifolds 

We start with a general remark. 

Remark 3.1. Assume that (^: M ^ M>o fulfills the assumptions of Proposition 2.1 Hi). 
Let be the sequenee produeed by Algorithm 1. Then we know that G 

(O, (^"(0+)/2]^^ where p — if v — \ (anisotropie ease), and p — )fQ, if v — 2 

(isotropie ease). By eonstruetion we have for the iterates produeed by Algorithm 1 that 

so that the sequence is monotonieally deereasing. By ( 8 ) and sinee (f is 

nonnegative, the funetion is bounded from below by zero and the sequenee {^Jiy{u^^\v 
eonverges to some b^. This holds also true for by (16). If M. is eompaet as in 

the ease of spheres or SO{3), then is elearly bounded. If AA is an Hadamard spaee 

as defined in the next subseetion and (p is eoereive, then is bounded sinee Jiy is by 

Proposition 3.3 eoereive. In these eases is also bounded and therefore there 

exists a subsequenee | }^. whieh eonverges to a point {u,v) G Mx [0,¥."(0)/2]^ 

For Hadamard spaces many results on the convergence of algorithms carry directly 
over from the Hilbert space setting. This is in particular true for the half-quadratic 
minimization algorithm. In this section we summarize these results for convex functions 
(p and data in Hadamard spaces. 

We start by recalling some basic facts. A curve 7 : [0,1] ^ X in a metric space (X, d) 
is called a geodesic if for all ti,t 2 G [0,1] the relation 

<^( 7 (^ 1 ), 7 (^ 2 )) =1^1 - t2M(7(0),7(l)) 

holds true. A function /i : X ^ M is called convex if /i o 7 is convex for each geodesic 
7 : [0,1] ^ X, i.e., if for all t G [0,1] we have 

h{^{t)) < i/i(7(0)) + (1 - ^)/^(7(l)) 

and strictly convex if we have a strict inequality for all 0 < t < 1. An Hadamard space is 
a complete metric space {%, d) with the property that any two points x, y are connected 
by a geodesic and the following condition holds true 

d{x, + d{y, < d{x, + d{y, + 2d{x, y)d{v, w), (23) 

for any x,y,v,w G X. Inequality (23) implies that Hadamard spaces have nonpositive 
curvature [3, 39] and Hadamard spaces are thus a natural generalization of complete 
simply connected Riemannian manifolds of nonpositive sectional curvature, the so-called 
Hadamard manifolds. For more details, the reader is referred to [5,26]. Unfortunately, the 
spheres and the rotation group are not Hadamard manifolds, while the symmetric positive 
definite matrices have this nice property. The following facts can be shown similarly as in 
M^, see [5, Lemma 2.2.9], [48]. 

Lemma 3.2. Let {PL, d) be an Hadamard spaee and F: PL ^ {+^} eonvex Ise 

funetion whieh is eoereive, i.e., satisfies F{x) -hoc whenever d{x,X{f) -hoc for some 
xo ^ PL. Then F has a minimizer. If F is eonvex, then any eritieal point is a global 
minimizer. If F is eoereive and strietly eonvex, then the minimizer is unique. 

In an Hadamard space {PL, d) we have that 



(Dl) d\ H X H ^ M>o and d?: % x% ^ M>o are convex, and 
(D 2 ) : 1-L M>o is strictly convex. 

Then we obtain the following proposition whose simple proof is added for convenience in 
Appendix A. 

Proposition 3.3. Let {M^d) — (Li^d) be an Hadamard manifold. 

i) Let (f: M>o ^ M>o and in the ease V ^ G, further assume that (f is eoereive. Then 
the funetions Jy, z/ = 1, 2, in (3) and (4) are eoereive so that they have a minimizer. 

ii) If in addition ip is inereasing and eonvex, then the funetions Jy, z/ = 1,2, are eonvex. 
If in addition V — G or cp is strietly eonvex, then the funetions Jy, z/ = 1 , 2 , are 
strietly eonvex and have unique minimizers. 

Under the assumptions of Proposition 2.1 hi), we have in our algorithm that > 0. 
Then, we see similarly as in the proof of Proposition 3.3 that the functionals 
V — 1,2, are coercive and strictly convex. Thus the minimizer exists and is unique. 

Theorem 3.4. Let {M.,d) = {TL,d) he an Hadamard manifold. Let ip:M. ^ M>o be an 
even, eontinuously differentiable, eonvex funetion whieh fulfills 

i) ifiVi), t>0 is eoneave, 

ii) limt^oo ^ ^ 0, 

iii) ¥?"(0+) > 0. 

In the ease V G we further assume that ip is strietly eonvex. Then the sequenee }ken 
generated by Algorithm 1 eonverges to the minimizer of Jy, z/ = 1 , 2 . 

The proof which follows standard arguments is given in the appendix A. Note that the 
assumptions of the theorem are fulfilled for the first two functions in Table 1. 

Remark 3.5. The assumptions in the theorem inelude those in Proposition 2.1. Note 
that the eonvexity of ip and (^"( 0 +) > 0 imply that ip'{t) > 0 for all t > 0. Additionally the 
eontinuity of ip' is required to make the funetion s in (11) eontinuous and the (striet) eon¬ 
vexity of ip to make the objeetive funetion strietly eonvex. Sinee ip is eonvex, its derivative 
is inereasing. Together with p)"{0-\-) > 0 this implies that ip'{t) > 0 for all t > 0 so that ip 
is inereasing for t > 0 and eoereive. 

4 Numerical Examples 

In this section we demonstrate the performance of Algorithm 1 for the functions ip from 
Table 1. These functions are known for their edge-preserving properties. Note that ipi was 
used in the “lagged diffusivity fixed point iteration” [53] for real-valued images and in the 
iteratively re-weighted least squares method [24] for S^-valued and P(3)-valued images. 
The function ip 2 is a Moreau envelope of the absolute value function, also known as the 
Huber function. Both functions are convex. The non-convex function ips was used for 
edge-preserving restoration of real-valued images in [16,33]. 

Unless stated otherwise we use the anisotropic approach and Newton’s method in our 
implementations. Although neither the spheres nor the rotation group are Hadamard 
manifolds we have observed convergence in all our numerical examples. This may be 
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ip{t) 

s{t) 

^l{t) 


1 

^2{t) 

{ t < e. 

f 1 t < £^ 

t t < £ 

1 # 

^3{t) 

1 — exp(—s^t^) 

exp(—s^t^) 


Table 1. Functions ip fulfilling the assumptions of Proposition 2.1. 


due to the fact that neighboring image pixels have values which are close enough on the 
manifold. 

The algorithms where implemented in MatLab Version 14b. The computations were 
performed on a Dell with 8 GB of RAM and an Intel Core i7, 2.93 GHz, on Ubuntu 14.04 
LTS. 

4.1 valued data 

We start with the one-dimensional signal in Fig. 1 to show how the different functions p 
from Table 1 perform and how the parameter s influences the results. The original signal 
in Fig. 1 (a) was obtained from f{x) — Sttx^ by sampling with size 0.01 and unwrapping 
modulo 27r such that the data are represented in [—7r,7r). Then wrapped Gaussian noise 
with standard deviation a — 0.3 was added. Using pi to restore the signal gives relatively 
larges error in Fig. 1 (b). The Huber function p2 and the exponential function p^ show 
better results in Fig. 1 (c) and 1 (d). The regularization parameter A was adapted to get 
the best error 

1 

err := — rf (/*,«*)> 

where N — 101. Making s in the Huber function larger leads to the smoother result in 
Fig.l (e) which approximates the original signal only well at the beginning of the signal. 
In Fig. 1 (f) we choose a larger 6: in the function p^ with the effect that edges of smaller 
height are smoothed and we have a staircasing effect for nearly equally high ascents. 

Next we want to demonstrate the difference between the anisotropic (12) and isotropic 
(13) half-quadratic minimization methods. To this end, the function atan2(x,^) was 
sampled over a regular grid with grid size resulting in Fig. 2 (a). Then we 

corrupt the image by removing a circular region from the center as shown in Fig. 2 (d). 
Using the anisotropic functional leads to Fig. 2 (b), where we observe artifacts in vertical 
and horizontal directions. The image produced by applying the isotropic functional in 
Fig. 2 (c) does not have this problem. This effect is also illustrated by the error plots in 
Figs. 2 (e) and 2 (f). 

4.2 valued data 

In our first example we denoise color images in the chromaticity and brightness space. For 

an RGB image the brightness is given by the real positive numbers b (i?^ + G‘^ + H^) ^ 
and the chromaticity by the S^-values c := (i?, G, B)/b. We want to mention that the first 
two examples consider problems which have values on the positive octant; only in the last 
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(a) Original and noisy signal 




(e) (/92, too smooth 



(f) staircasing 


Figure 1. Restoration of a noisy cyclic signal by half-quadratic minimization with various 
functions (p. (a) original (red) and noisy (blue) signals. Restored signal (blue) using (b) 
pi, 5 = 6 X 10“^, A = 3.4, err = 0.1007, (c) p2, 5 = 5 x 10“^, A = 5.2, err = 0.1007, 
(d) (^ 3 , 5 = ^, A = 10, err = 0.1001, (e) p 2 , e = I, \ = 20, err = 0.1733, (f) 

(/93, 5 = \/5, A = 10, err = 0.3756. 
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(a) Original image. (b) Anisotropic model. (c) Isotropic model. 



(d) Corrupted image. (e) Error of (b). (f) Error of (c). 


Figure 2. Inpainting of an image with cyclic data using the anisotropic and the isotropic 
model, (a) original image, (d) corrupted image. Restoration with (b) the anisotropic 
model and (c) the isotropic model using the function with A = 0.001 and 5 = 10“^. 
(e) and (f) error images. 
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(a) Original image. 


(b) Noisy, PSNR: 20.31. 


(c) ifi, PSNR: 30.29. 



Figure 3. Denoising in the chromaticity-brightness space, (a) Original image “Peppers”, (b) 
Corrupted image by Gaussian noise on RGB, a = 0.1. Restored images by (c) ifi using 
5c = 10“^, Sb = 10“^, Ac = 0.44, = 0.08, (d) ip 2 using 5c = 10“^, 5 ^ = 10“^, Ac = 

15, A^ = 10, (e) (ps using 5c = 2 V^, 55 = \/^, Ac = 0.1, A 5 = 0.03, and (f) TV method 
in [56], a = 0.05. 


example we look at data covering the whole sphere. We compare half-quadratic minimiza¬ 
tion with the different functions p and the TV approach from [56] in Fig. 3. We took the 
image “Peppers”^ in Fig. 3 (a) and added Gaussian noise with standard deviation a = 0.1 
to all three color channels in the RGB model. For denoising in the chromaticity-brightness 
model we optimized A with respect to best PSNR for both channels separately using a grid 
search on for pi and ps^ and for p 2 on N. Furthermore, we optimized 5 first in order 
of magnitude k = 10-^ and refined the search on ^N. Both channels are restored using 
the same function. The square-root functional pi shows smoother transitions at edges 
while the Huber function p 2 tends more to staircasing. Using the exponential function p^ 
yields a worse PSNR and hence does not compete with the previous two functions; the 
edges are preserved but within the more constant regions some noise is left. The bright 
spot (see upper magnification) appears also too smooth. This originates from the too 
smooth transitions in the brightness which are not detected as edges. The TV regulariza¬ 
tion introduces staircasing and it is not able to reduce the noise in the dark area (lower 
magnification). 

In the second example we use half-quadratic minimization for colorization in the 
chromaticity-brightness space. We assume that the brightness of the image is known, 
but 99 percent of the chromaticity information is lost. The original image is shown in 
Fig. 4 (a) and its corrupted version in Fig. 4 (b). For inpainting the chromaticity we have 
used a nearest neighbor initialization. With the regularizing function pi we obtain the 
result depicted in Fig. 4(c). We compare this with Fig. 4(d) which is obtained by using 
the chromaticity colorization method in [37] which we have implemented for comparison. 

Our final experiment shows the smoothing of 3D directions in the synthetic image in 

^Taken from the USC-SIPI Image Database, available online at http://sipi.use.edu/database/ 
database. php?voluine=misc&image=15 
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(a) Original image (b) 99% color lost. (c) PSNR: 27.19. (d) [37], PSNR: 22.49. 

house. 


Figure 4. Image colorization. (a) Original image, (b) Corrupted image where 99% of 
the color information (chromaticity) is lost. Colorization using (c) inpainting of the 
chromaticity with A = 1, 5 = 10“^. (d) the method in [37] with parameters 

r = l,p = l,cri = 2,(72 = 00,7 = 0. 



(a) Original S^-field. (b) Noisy S^-field. (c) (^ 2 - 


Figure 5. (a) Original S^-field of size 64x64. (b) Corrupted field by Gaussian noise, a = 0.1. 

(c) Restored field with ip 2 , X = 2.6, 5 = 10“^, leaving an error err = 0.1705. 

Fig. 5 (b). We use half-quadratic minimization with (p2 to obtain Fig. 5 (c). The original 
pattern is again visible. 

4.3 P(3)-valued data 

Our first example illustrates the inpainting capabilities of the half-quadratic minimization 
method by an artificial example. The 'P(3)-valued image of size 16 x 16 in Fig. 6 (a) has 
a jump at I in X direction. We destroy a center square of size 12 x 12, see Fig. 6 (b). We 
are able to reconstruct the inpainting area nearly perfectly by using cpi and 6: = 10“^, see 
Fig. 6(c). Nevertheless decreasing 6: introduces more and more staircasing, cf. Fig. 6(d). 
This resembles the TV case shown in Fig. 6 (e), using the model from [56]. 

An important application of 7^(3)-valued image denoising is Diffusion Tensor Magnetic 
Resonance Imaging (DT-MRI). The Camino project^ [17] provides a DT-MRI dataset of 
the human head which is freely available.^ The complete data is given as a 3D image 
/ = {fij,k) C 7:>(3)112x112x 50^ where we apply the half-quadratic minimization to each 
of the traversal planes k G {1,...,50}. The original dataset, cf. Fig. 7(a), is plotted 

^see http://cmic.cs.ucl.ac.uk/camino 

^follow the tutorial at http://cmic.cs.ucl.ac.uk/camino//index.php?n=Tutorials.DTI 
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(a) Original image. 



(b) Corrupted image. 



Figure 6. Inpainting of an 7^(3)-valued image, (a) Original image, (b) Image with unknown 
areas. Half-quadratic based inpainting (c) yields a perfect result. Decreasing 5 like in (d) 
yields a result closer to (e) TV. 



(a) Original “Camino” data set. (b) a = 0.1, 5 = 10 

Figure 7. Denoising with half-quadratic minimization of every traversal plane of the 
“Camino” data set (DT-MRI of the human head). 
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(d) Grain. 


Figure 8. (a) The raw EBSD data of a Magnesium sample, (b) the colorization of the sphere 
used to assign to each rotation a certain color according to the mapping SO(3)/5' 3 m 

1)^ G S^/5', (c) colorization of a spherical triangle, (d) stretched colorization 
for one grain. 


using the anisotropy index relative to the Riemannian distance [32] normalized onto [0,1) 
and colored in hue. The half-quadratic minimization is used with (pi and the parameters 
A = 0.1, s = 10“^ and a maximum change between two successive iterations being 10“^^ 
as a stopping criterion. We obtain the result shown in Fig. 7 (b). For the complete dataset 
of 168,169 nonzero matrices, the algorithm needed 2,492 seconds to compute the result. 

4.4 SO(3)-valued data 

Processing images with SO(3)-valued entries is fundamental in the analysis of polycrys¬ 
talline materials by means of Electron Backscattered Diffraction (EBSD), cf. [2,28]. Since 
the microscopic grain structure affects macroscopic attributes of materials such as ductil¬ 
ity, electrical and lifetime properties, there is a growing interest in the grain structure of 
crystalline materials such as metals and minerals. EBSD provides us for each position on 
the surface of a specimen with a so called Kikuchi pattern, which allows the identification 
of the structure (material index) and the orientation of the crystal at this position relative 
to a fixed coordinate system (SO(3) value). Since the atomic structure of a crystal is 
invariant under its specific symmetry group S C SO(3) the orientation is only given as an 
equivalence class [mo] = {mos | s G 5} G SO(3)/5, mo G SO(3). 

Fig. 8 (a) displays a typical EBSD image consisting of lattice orientations of deformed 
Magnesium collected by [43]. Each pixel of the image corresponds to a position on the 
surface of a Magnesium specimen. The color of the pixels is chosen corresponding to the 
orientation measured at this position according to the following color mapping: for a fixed 
vector r G we consider the mapping $: SO{3)/S -3^ [m] [m~^f]. Next we 

colorize the quotient as it is depicted in Fig. 8(b). From the colorization scheme 

the symmetry group S C SO(3) of Magnesium becomes visible which has six rotations 
with respect to a 6-folded axis (A:7r/3, k = 1,... ,6 rotations around c direction) and six 
rotations with respect to 2-folded axis ai, a 2 perpendicular to that. 

EBSD images usually consist of regions with similar orientations called grains. For 
certain macroscopic properties the pattern of orientations within single grains is of im¬ 
portance [7], e.g., for the computation of geometrically necessary dislocations [35,49] the 
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gradient of the rotations within single grains has to be determined. As the rotation deter¬ 
mination by Kikuchi patterns is sometimes fragile, the rotation valued images determined 
by EBSD are often corrupted by noise and suffer from missing data so that denoising and 
inpainting techniques have to be applied [25]. For detecting grains in the raw EBSD data 
we applied a thresholding algorithm [ 8 ]. Fig. 9 (a) displays a single grain with its rotations. 
Since the rotations vary very little within a single grain we applied a sharper colorization, 
cf. Fig. 8 (d) to make the noise and the rotation gradient visible. 

We want to apply half-quadratic minimization to denoise EBSD images. Since crystal¬ 
lographic symmetry groups are finite the quotient SO{3)/S is locally isomorphic to SO(3). 
In particular, the formulas given in Appendix B.3 can be applied. In a first experiment 
we apply half-quadratic minimization using (fi to the rotation-valued image depicted in 
Fig. 9 (a) which leads to the smooth image in Fig. 9 (b). In a second experiment we ran¬ 
domly removed 30% of the data shown in Fig. 9 (c). Using half-quadratic minimization for 
jointly inpainting and denoising the image we obtain the result shown in Fig. 9 (d) which 
looks very similar to those in Fig. 9 (b). In a third experiment we applied half-quadratic 
minimization simultaneously to several grains. The challenge from the mathematical point 
of view is that SO(3) is not an Hadamard manifold. Convergence can be guaranteed only 
locally, which is the case for single grains but may be not true when considering several 
grains simultaneously. From the practical point of view this case is especially interesting 
as missing data usually occur at grain boundaries, i.e., between grains. However, for our 
data we have got promising results. Fig. 10 (a) shows the grain from the previous example 
(pink color) and two other grains in its neighborhood. Note that the top middle area (light 
green) and top right area (brown-green) belong to the same grain. Pixels with missing 
data are plotted white. Half-quadratic minimization restoration with (pi improves the im¬ 
age as can be seen in Fig. 10 (b). For our last experiment we again randomly remove 30% 
of the data, cf. Fig 10(c). Restoration with ipi leads to the result in Fig. 10(d), which 
is again hardly to distinguish from Fig. 10(b). Using (^3 leads to even better results as 
depicted in Fig. 10 (e). We can adjust the smoothing in such a way that the edge distin¬ 
guishing two grains is not smoothed, while smaller rotation changes are smoothed. This 
is advantageous in the large top grain. Finally Fig. 10 (f) shows the zoom to the (pink) 
grain from Fig. 9 with adapted color map. 

5 Conclusions 

We adapted the principles of half-quadratic minimization to the setting of complete, con¬ 
nected Riemannian manifolds. In particular, the notation of the c-transform provides an 
interesting point of view. For Hadamard manifolds we proved the existence and unique¬ 
ness of the minimizer of the corresponding functionals as well as the convergence for the 
alternating minimization algorithm under moderate assumptions. The multiplicative half¬ 
quadratic minimization method resembles a quasi-Newton method [33] and appears to be 
very efficient in our numerical examples. There are numerous applications of the approach. 
In this paper, images having values in a manifold such as the 2-sphere or the symmetric 
positive definite matrices were denoised. The method was also used for inpainting missing 
information into images consisting of either rotation matrices or symmetric positive def¬ 
inite matrices. In the chromaticity-brightness color model, the inpainting technique was 
applied to the task of colorization. The method has further potential in EBSD. 

Topics of future research are the derivation of convergence proofs for more general 
manifolds under special assumptions on the local behavior of the data. Furthermore, 
different data terms must be included for other applications, and the inclusion of higher 
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(a) Original grain. 


(b) Grain smoothed with (fi. 



(c) Grain with 30% lost data. 



(d) Grain restored with ifi. 


Figure 9. (a) Original grain from a Magnesium specimen (grain boundary in black), (b) 
smoothed grain using ifi with parameters A = 0.1, 5 = 10“^, (c) grain with 30% lost data, 
marked in white, (d) inpainted and smoothed grain with ifi and the same parameters. 


order differences into the regularization term of the model is of interest. 

A Proofs 

Proof of Proposition 2.1. 1. In the additive case we have 

=si - f(')}+ 

= — sup| st — — (fit)) I + 

*(=R '--- '' 2a 


1 


By assumption on $ we know that $ = which implies 

= inf {-st + = inf {-ts + $*(«)} + 

= —= ^{t)- 

This finishes the proof of i). The function 

h{t) := c{t^ s) — (p(t) = -at^ — (p{t) — st + — 5 ^ 


18 


























(a) Three grains. (b) Grains smoothed with ipi. (c) Grains with 30% lost data. 



(f) Grain from Fig. 9. 


Figure 10. (a) Three grains from a Magnesium specimen with missing data in white, (b) 
Smoothed (and ipainted) grains with ipi and parameters A = 0.15 5 = 10 “^. (c) Grains 
with 30% data lost. Restored grains from (c) with (d) ipi and parameters A = 0.05 
5 = 0.5 X 10“^, (e) (^3 and parameters A = 0.1, 5 = v^SO- (f) Grain from Fig. 9 after 
restoration with 993 (original boundary in black). 


is continuous, convex and by (10) coercive so that the global minimizer in (7) is attained 
for 0 = h'{t) — at — ^'{t) — 5, i.e., for (t, at — ^'{t)) which proves ii). 

2. In the multiplicative case we obtain, since (f is even, 


■= inf = inf{f^s - 


t>0 ' 


= inf{rs-v?(x/r)} = 


sup 

r>0 


|-rs - 


We have 


SO that we can restrict our attention to t > 0. By assumption, $ is convex and Isc. Thus, 
$ = $** and we obtain for t > 0 that 


— inf {ts — (^^"^( 5 )} = — sup{—— $*(—s)} = —$(t) = (^(v^). 


This yields i). To see ii) we first note that condition (9) implies for 5 > 0 that the 
objective function in (7) is coercive such that the infimum is attained. For s < 0, we have 
(■p^{s) = —oc. For 5 > 0, we obtain 


arg minlt^s — 
t>o 


V’W} = ( 


arg minlrs — ip 

r>0 



By assumption on the function h{r) := rs — ip{^) is convex in M>o. A global minimizer 
of h is attained either for the solution of 0 = h'{r) — s—i^ip'{^) if this solution is positive 
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or for r = 0 if lim^^o+^ 0, i.e., 5 > ^(^"(0+). Therefore (t, 5 ) = (O, ^(/^"(0+)) is a 
solution. 

Finally, the concavity of for t > 0 implies that (p'{\/i)/{2\/i) and thus (p{t)/{2t) 

is decreasing. Under the additional assumption in iii) we get s G (O, ^ 2 ^^^ ] • D 

Proof of Proposition 3.3. i) If V = let /)||2 ^ 00 , where d{u^ f) := (^d{ui^ fi)).^g. 
Then Fiu) i d^{ui^fi) goes to infinity and the functionals J^y, z/ = 1,2, are 
coercive. In the case V ^ Q choose io ^ V and let izq be the constant image with entries fi^. 
Let ||d(rz,izo)||2 ^ 00 . Assume that Jjy{u) remains finite, so that in particular (i(/^Q,iz^Q) 
and d{ui^ Uj)^ j G A/’(i)^, i ^ G are finite. By the construction of the neighborhoods J\f{i)^ 
there exists for every j ^ G ^ path io, A,..., = j with ii^i G J\f{ii)^ and 

d{fi^,Uj) < d{fi^,Ui^) + d{ui^,Ui^) + .. . + d{uij^_^,Uj). 

Since the right-hand side remains finite this contradicts ||d(rz, izq)II 2 ^ 00 . Hence 
z/ = 1, 2 are coercive. 

ii) By (D2) we have that F{u) is convex and strictly convex if V = 0. If a function 
h: ^ M is convex, then, for any geodesic 7 : [0,1] ^ joining G we obtain 

since (p is increasing and convex that 

p) o h(pf{t)) < p(th{x) + (1 — t)h{y)) < t{p o h){x) + (1 — t){p o h){y) (24) 

so that p o h is convex. If p is strictly convex the last inequality is strong so that p o h 
is strictly convex. With h := hij = d{ui^Uj) : Td? M this implies by (Dl) that Ji 
is convex, resp. strictly convex. Concerning J 2 notice that the convexity of /ii(xo,x^), 

i = 1 ,..., — 1 , on Td? implies convexity of /i(xo,..., h,^-i) := ^ on 

by 

K, — l K, — l 2 

(7o(i),7i(i)) < ^(i^i(7o(0),7i(0)) + (1 - t)/ii(7o(l),7i(l))) 

i=l i=l 

K—1 

= rh^ijoiO)) + {I - {joir) +2t(t- l)y]/ii(7o(0),7i(0))/i,i(7o(l),7i(l)) 

i=l 

and by the Schwarz inequality 

< rh^{^o{0)) + (1 - t)^/i^( 7 o(l)) +2t{t- l)/i( 7 o( 0 ))/i( 7 o(l)) 

= (t/i,(7o(0)) + (1 - t)/i(7o(l))) . 

In the case V — G^ strict convexity follows by the strict convexity of the data term 
^ieQ strictly convex p by the strict convexity in (24). □ 

Proof of Theorem 3.^. By Remark 3.1 we know that lim^c^oo =: b and that 

there exists a subsequence |which converges to some {u^v). Since is 
continuous we have lim^^oo = Jiy{h^v) = b. Let 

V := s (du) = argmin J7^(iz,'z;) and u := argmin^jy(iz, i)). 

V u 

The continuity of 5 and d implies that lim^^oo = lim^^oo '^(^^( 70 ) ~ ^{^u) = ^ and 

the continuity of T{u) := argmin^ J7’iy(2;, 5(d^i)) that lim^-^oo= li^ij^oo = 

T{u) = u. By (22) we conclude 

b = lim = Ju{u,v) < Ju{u,v) < X{u,v) = h. 

j^oo 
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Thus, Jiy{u^v) = J^{u^v) — J^{u^v) and since and have unique min- 

imizers we obtain that u = u and v = v. Consequently, v = s{du) and u — u — 
argmin^ J^{u^ s{du))^ which by Remark 2.2 implies that iZ is a minimizer of Jjy. 

Assume that lim^^^oo A Then there exists 6 : > 0 such that infinitely many 
not contained in the open ball B^{u) of radius s centered at u. Since is bounded, 

there exists a convergent subsequence which converges to some point xx* ^ xx in the 

closed set M\Bs. Then limj^oo Ji/(xx^^"-^^) = J^iu") — J{u) which contradicts the fact that 
Jy has a unique minimizer. □ 

B Exponential and Logarithmic Maps 

B.l The Sphere 

We use the parametrization 


x{e,ip) 


^cos if cos 9^ 
sin ip cos 9 
sin 0 


[ 0 , 27 r). 


Then we have the tangent spaces 

r^(§ 2 ) = Ta.( 0 _^)(§^) := [t] G : rpx = O} = span{ei( 6 »,y?),e 2 ( 6 »,y?)} 
with the normed orthogonal vectors 




cos if sin 9^ 
- sin (p sin 9 
cos 9 


62(9, (f) := 


1 dx 
cos 9 dp) 



The Riemannian metric is just the Euclidean distance in M^. The geodesic distance is 
given by d^ 2 {xi^X 2 ) := arccos(xi, X 2 ), and the exponential map and (locally) its inverse, 
resp., by 

71 

exp^{tr]) := cos(t||x7||2)x + sin(t||x7||2) 

1 ^2 (^X\^ X2')X\ . . 

log^ X2 := -n-^—rp arccos(xi,X2). 

\\X 2 — (3:1, X 2 )xi ||2 


B.2 The Manifold V{r) of Symmetric Positive Definite Matrices 

By Exp and Log we denote the matrix exponential and logarithm defined by 

00 ^ 00 ^ 

Expx:='^—x^, Logx - x)^, p{I - x) < 1. 

k=0 k=l 

Let Sym(r) denote the space of symmetric r x r matrices with (Frobenius) inner product 
and norm 

(25) 

Let V{r) be the manifold of symmetric positive definite rxr matrices. It has the dimension 
dimP(r) = n = tangent space of V{r) at x G V{r) is given by TxV{r) — 


{A^B) ciijbij^ \\A\\ i a^j 
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{x} X Sym(r) := {x 2 7 yx 2 : ry G Sym(r)}, in particular TjV{r) — Sym(r), where / denotes 
the r X r identity matrix. The Riemannian metric on TxV reads 

{m,'n2)x ■= tT{r)ix~^r]2X~^) = {x~^r]ix~^,x~^r]2X~^), 

where (*,•) denotes the matrix inner product (25). The geodesic distance is given by 
_ 1 _ 1 

d'p{xi^X 2 ) := ||Log(x]^ ‘^X 2 X^ ^)||, and the exponential map and its inverse by 

1 111 1 _1 _1 1 
exp^{tri) := x 2 Exp(tx~^r]x~^)x^ ^ resp., log^,^ X 2 := x^ Log(^x-^ ‘^X2X^ 5 

see [44]. 

B.3 The Manifold SO(3) of Rotation Matrices in 

The manifold of 3 x 3 rotation matrices is defined as 

SO(3) := {x G I x^x = I and detx = l}. 

The geodesic distance between two rotation matrices xi^X2 G SO(3) is given by 

<^SO( 3 )(^i 5 ^ 2 ) — \/2arccos 

The tangential space at x G SO(3) reads 

r^SO(3) := {xv:ve r/SO(3)}, r/SO(3) := {v G :v + v^ ^O}. 

For 7] G Ta;SO(3) the exponential map at x G SO(3) and (locally) its inverse are defined 
as 

exp^(77) := xExp(x^77), resp., log^,^ X 2 = xiLog(x^X 2 ). 

The SO(3) can be parametrized in various ways. Due to the form of our data we prefer 
to use quaternions for the representation and the similarity of SO(3) to the group see, 
e.g., [23]: We decompose the unit quaternions q — ( 5 ,'?;^)^ G into a real part 5 G M 
and a vector part 'y G M^. The multiplication of two quaternions ^ 1,^2 G is given by 

f Sl\ f S 2 \ f SiS 2 - V1V2 

qi o q 2 := { o = 

\vi J \V 2 J \S2V1 + S2V1 +vi X V 2 

with the conjugate quaternion q := (5, —as inverse element and unit element (1,0,0, 0)^ 
The quaternions can be identified with the rotations of SO(3), where the quaternions q 
and —q correspond to the same rotation. More precisely, is a double cover of SO(3), 
see, [12, Chap. Ill, Sect. 10]: 

We work with the representative having a positive first component. With this representa¬ 
tion at hand, he Riemannian metric of SO(3) can be deduced from the Euclidean metric 
in M^. The geodesic distance can be written as 

d^3{qi,q2) = 2arccos |(gi,g 2 )|. 



^ 1 - tr(x^X2) ^ 
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The exponential map of 77 G TpS>l is given by 


exp T] '.= sgn s 


, := gcos||r ?||2 + 


■ Sin 


and the logarithmic map at qi of q 2 by 

92 - (92,91)91 


loggi 92 := 


I92 - (92,91)91 II2 


arccos(|(gi,g 2 )|) sgn((gi,^ 2 ))- 
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